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ABSTRACT 

We investigate the influence of the cooHng epoch on the formation of galaxies in 
a cold dark matter dominated universe. Isolated haloes, with circular speeds typical 
of spiral galaxies, have been selected from a low resolution numerical simulation for 
re-simulation at higher resolution. Initial conditions for each halo consist of a high 
resolution region containing dark matter and gas, and a nested hierarchy of particles 
representing the mass distribution of the original low resolution simulation. These 
initial conditions are evolved with two smoothed particle hydrodynamics (SPH) codes, 
TREESPH and GRAPESPH, so that discrepancies due to differences in evolutionary 
and star formation algorithms can be analysed. In previous SPH simulations, strong 
outward transport of angular momentum has led to the formation of disc-like systems 
with much smaller angular momenta than observed in real disc galaxies. Here we 
investigate whether this problem can be circumvented if feedback processes prevent 
disc formation until late epochs. In some of our models, the gas is evolved adiabatically 
until a specified redshift Zcooi, at which point the gas is allowed to cool radiatively 
and star formation may begin. In other models, cooling is continuous throughout the 
evolution but suppressed by a factor chosen to allow a discs to grow roughly linearly 
with time. The results of varying the cooling epoch for each of five different haloes 
are analysed. When cooling and star formation are initiated at redshift Zcooi — 4, 
stellar discs are destroyed during merger events and we observe similar catastrophic 
transport of angular momentum as seen in previous work. With cooling suppressed 
until Zcooi = 1 , discs can form by the present day with angular momenta comparable to 
those of observed disc galaxies. We conclude that feedback processes, which prevent gas 
from collapsing until late epochs, are an essential ingredient in disc galaxy formation. 



1 INTRODUCTION 



Cosmological N-body simulations have shown that the evo- 
lution of small density perturbations in an expanding uni- 
verse can lead to a distribution of matter very similar to 
the observed large-scale structure (see e.g. Jenkins et al. 
1998). Purely N-body simulations have also led to an un- 
derstanding of the evolution of dark matter on the scale of 
individual galaxies. For example, simulations of cold dark 
matter (CDM) models have shown how dark matter haloes 
form by the hierarchical merging of sub-units (Frenk et al. 
1988), acquire angular momentum via tidal torques (Barnes 
and Efstathiou 1987), and develop radially-averaged density 
profiles that satisfy simple scaling relations (Navarro, Frenk 
& White 1997). 

Gas-dynamical simulations of the formation of indi- 
vidual galaxies within dark matter haloes have been much 
less successful. In simulations without star formation {e.g. 
Navarro & Steinmetz 1997), the gas component can shock 
and dissipate energy (unlike the coUisionless halo material) 
forming thin discs within the dominant halo sub-clumps. 
The gas discs may be disrupted in subsequent mergers, but 



the efficient conversion of kinetic energy to thermal energy 
and rapid cooling allows the gas to settle into a new disc on 
a timescale that is short compared to the Hubble time. The 
evolution thus proceeds through a hierarchy of disc forma- 
tion and merging, leading, in the absence of star formation, 
to gaseous discs at the present day. However, during merg- 
ing, angular momentum is efficiently transported outward 
into the halo and the resulting gaseous discs have angu- 
lar momenta some two orders of magnitude below those of 
observed spiral galaxies (Navarro & Steinmetz 1997). Evi- 
dently, such simulations cannot account for the formation of 
real disc galaxies. We will sometimes refer to this problem 
as the 'angular momentum catastrophe'. 

On the other hand, in simulations in which efficient star 
formation is included, stars form within dark matter sub- 
clumps at early times. These largely stellar sub-units then 
merge as the system evolves, losing angular momentum and 
forming a hot spheroidal component by the present day. The 
final objects in such simulations are elliptical-like with no 
strong disc structure {e.g. Katz 1992). 

The failure to produce realistic disc galaxies in numer- 
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ical simulations contrasts with simple analytical models of 
disc galaxy formation. These models, first developed by Fall 
and Efstathiou (1980) and Gunn (1982), show that the sizes 
and collapse times of disc galaxies can be understood wi- 
hin the context of a two-component (gas and dark matter) 
theory of galaxy formation (White and Rees 1978) provided 
the gas component conserves its angular momentum during 
collapse. More recent applications of this type of model have 
been developed to explain other properties of disc galaxies, 
e.g. the slope and scatter of the TuUy-Fisher (1977) relation 
(Mo, Mao & White 1998) and the distribution of disc sur- 
face brightnesses (Dalcanton, Spergel and Summers 1997). 
In fact, we will show in Section 2 of this paper, that simple 
arguments can be developed to show that in CDM mod- 
els, galactic discs cannot have lost much angular momen- 
tum during collapse. There is, therefore, an apparent incon- 
sistency with the 'angular momentum castrophe' described 
above. 

In this paper we suggest that this conflict might be 
solved if stellar feedback processes prevent discs forming un- 
til late epochs [z ^ 1-2). If some fraction of the protogalac- 
tic gas is prevented from collapsing until late epochs, when 
galactic dark matter haloes are evolving slowly and most of 
their substructure has been erased, then it may be possi- 
ble to form discs without the angular momentum transport 
characteristic of merging sub-units. Our aim in this paper is 
to test this idea using SPH simulations. 

We review the observed angular momenta of disc galax- 
ies in Section 2 and compare them with the angular mo- 
menta of CDM haloes determined from N-body simulations. 
We discuss also the surface brightness distributions of disc 
galaxies which we relate to the large dispersion in the angu- 
lar momenta of CDM haloes. The SPH codes that we use, 
and various other numerical details such as the star forma- 
tion algorithms, are described in Section 3. The generation 
of initial conditions for the SPH simulations is described in 
Section 4 and our main results are presented in Section 5. 



2 THE ANGULAR MOMENTA OF SPIRAL 
DISCS 

In this section, we review observational constraints on the 
angular momenta of disc galaxies and compare the results 
with the angular momenta of dark haloes in N-body simu- 
lations of a CDM model. 

We first investigate the angular momenta of haloes iden- 
tified in numerical simulations of aii fl — 1 scale invariant 
CDM model. There have been many previous investigations 
of this problem, including Barnes and Efstathiou (1987), 
Frenk et al. (1988), Zurek, Quinn and Salmon (1988), and 
more recently by Lemson and Kauffmann (1997). The re- 
sults described here are based on the simulation that is used 
to generate the initial conditions of the high resolution SPH 
simulations described in the main part of this paper. The 
numerical simulation used here is run (d) of the ensemble de- 
scribed by Efstathiou (1995). The simulation is of a spatially 
flat SI = 1 CDM model with scale-free Gaussian initial condi- 
tions generated from the power spectrum given by equation 



(7) of Efstathiou, Bond & White (1992) with F = flh = 0.^ 
We assume h = 0.5 in the rest of this paper, unless explicitly 
stated otherwise. The simulation contains A'p = 10® parti- 
cles within a cubical volume of box-length Lbox = 40 Mpc, 
with the spectrum normalized so that the rms mass fiuc- 
tuations in spheres of radius 8/i~^Mpc is ag = 0.59 at the 
present day. This normalization approximately matches that 
inferred from the abundances of rich clusters of galaxies in 
an Q, — 1 universe (White, Efstathiou & Frenk 1993; Eke, 
Cole & Frenk 1996) . The mass per particle in the simulation 
is 4.43 X 10^ Mq. 

Haloes are identified by linking together particles sep- 
arated by 0.2 times the mean interparticle separation using 
the percolation algorithm described by Davis et al. (1985). 
The mean overdensity of these haloes is A « 200, corre- 
sponding roughly to the critical overdensity defining the viri- 
alized region of a collapsed uniform sphere. For each halo, 
we compute the angular momentum J, mass M, total energy 
E and circular speed Vc defined by 



= G 



M(< Rc 
Rc 



(1) 



where M(< Rc) is the mass contained within a sphere of 
radius Rc centred at the centre of mass of the halo. We 
used a radius of Rc — 0.4Mpc, which is approximately the 
virial radius of a halo with a circular speed of 200 kms~^ 
(see equation We also compute the dimensionless spin 
parameter A (Peebles 1969) 



-5/2 



G~ 



(2) 



Fig. |l^ shows the specific angular momenta for haloes 
identified at z = plotted against their circular speeds. 
This type of diagram is more instructive than the usual plot 
showing specific angular momentum against either halo or 
disc mass (Fall 1983, Navarro & Steinmetz 1997), because 
the circular speed of a collapsed, centrifugally supported disc 
is likely to be closely similar (to within 50%) of the halo 
circular speeds (see Section ^ and Mo et al. 1998). With 
a diagram such as Fig. ^, the circular speeds measured for 
real disc systems can be used as indicators of the circular 
speeds, and typical specific angular momenta, of the haloes 
in which they formed. 

For a halo with density profile 



(3) 



-3/2 



the virial radius is given by 

= 0.3(1 + z)''/'i;3oo(A2oo)"''''fe"'Mpc, (4) 

where we have assumed Q = 1 and the notation vsoo ~ 
i;c/300kms~i , A200 = A/200, etc. If we assume that the 
halo density profile truncates at r = Rh, we can compute 
the specific angular momentum from the spin parameter A 

j{RH) = V2\VcRh 

^ 6.5 X 10^(1 + zy^^^Xo.osVaooh''^ kms"^kpc"\ (5) 



* Where h is the Hubble constant in units of 100 kms ^Mpc ^ 
and Q is the cosmological density parameter. 
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Figure 1. (a) The circles show the specific angular momenta of haloes with an overdensity A si 200 identified in an f2 = 1 N-body 
simulation of a CDM universe plotted against their circular speeds. The solid line shows the specific angular momentum expected for a 
halo with dimensionless spin parameter A = 0.05 (equation (b) as (a) except that the specific angular momenta have been computed 
for haloes identified at z = 1 and the circular speeds are those of the haloes at z = within which they merge. The solid line is equation ^ 
with 2 = 1; (c) the symbols show the specific angular momenta of a sample of disc galaxies described in the text and the solid line shows 
equation H as plotted in panel (a). 



where we have chosen A — 0.05, since this is approximately 
the median value of A measured for haloes in CDM models 
(Barnes & Efstathiou 1987). Equation ^ is shown as the solid 
line in Fig. |l| and provides an excellent approximation to the 
median specific angular momentum as a function of circular 
speed. 

The results plotted in Fig. |l|a probably provide a conser- 
vative upper limit to the angular momenta of disc galaxies. 
The angular momentum of the disc component will depend 
on the epoch at which the gas that forms the disc can cool, 
and on its radial extent. The specific angular momentum 
of discs will be lower than shown in Fig. if either the 
disc gas collapses from well within the virial radius Rh, or 
if the disc collapses at high redshift. For example, Fig. |l|b 
shows the specific angular momenta of haloes identified at 
z = 1 plotted against the circular speed of the halo at z = 
within which they merge. The solid line in the plot shows 
equation |5| The specific angular momenta of haloes at z = 1 
are typically three times smaller than at 2 = 0, because the 
virial radius at 2 = 1 is about three times smaller than at 
z = and haloes typically display a nearly constant rotation 
speed with radius (Frenk et al. 1988). 

The third panel in Fig. |l| shows the specific angular mo- 
menta for a large sample of disc galaxies plotted against their 
circular speeds. The photometric parameters are from Byun 
& Freeman (unpublished). These authors have analysed I- 
band CCD images of galaxies from the Mathewson, Ford 
& Buchhorn (1992) sample, using a two-dimensional pro- 
file fitting algorithm described by Byun & Freeman (1995). 
Each image has been decomposed into r'^''*-law bulge and 
exponential disc components, providing estimates of the ex- 
ponential disc scale length , central surface brightness Iq, 
bulge-to-disc ratio B/D, inclination angle and I-band appar- 
ent magnitude. We use estimates of the circular speeds Vc 
given in Mathewson et al. (1992), corrected for inclination. 



to compute the specific angular momenta for a centrifugally 
supported exponential disc with a fiat rotation curve. 

^ = 2v.a-\ (6) 

We assume Hq = 50 kms~^Mpc~^ to compute distances and 
include only gala:xies with distances > 10 Mpc and bulge- 
to-disc ratios B/D < 0.3, to reduce errors arising from inac- 
curate distances and poorly determined disc profile param- 
eters. These selection criteria reduce the number of galaxies 
from 1355 in Byun & Freeman's sample to 1036. The solid 
line in Fig. file shows equation I with z = {i.e. as plotted 
in Fig. ^). This line lies about three times higher than the 
median defined by the data points. This is encouraging be- 
cause it shows that the angular momentum of a typical disc 
galaxy is similar to the angular momentum within the viri- 
alized region of a present day CDM halo of similar circular 
speed. 

Fig. |l| also shows that the scatter in the specific angular 
momenta of real disc galaxies is comparable to the scatter of 
the specific angular momenta measured for N-body haloes. 
This can be demonstrated using the following simple argu- 
ment (Efstathiou & Barnes 1984). Assume that a fraction 
fo of the baryonic mass within the virial radius of a halo 
collapses to form a disc. The luminosity of the disc is 

-^/-'"(ff)„^/"4(5)„a)"'-- 

where the expression for Mh has been substituted assum- 
ing that Q, = 1 and z = 0. The proportionality expresses the 
tight TuUy-Fisher (1977) correlation between luminosity and 
circular speed observed for real disc systems. Evidently, the 
(uncertain) cooling and feedback processes which determine 
fo must be tightly correlated with circular speed to repro- 
duce the TuUy-Fisher relation. The correlation between disc 
luminosity and circular speed for the Byun-Freeman sample 
is shown in Fig. ^a. A least squares fit to the data points 
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Figure 2. (a) shows the disc luminosity of the Byun &; Freeman sample plotted against circular speed i^c- The line shows a least squares 
power law fit Lc oc with a = 2.67. (b) compares the distribution of lov^ for the Byun & Freeman sample (solid histogram) with 
the distribution of determined for haloes in the CDM simulation, scaled to have approximately the same mean as the observations, 

(dotted histogram) and the best-fitting lognormal distribution for A (solid line). If the disc component conserves angular momentum 
during collapse then Iqv^ should be distributed as 



gives an exponent a = 2.67 ± 0.04. If we equate the spe- 
cific angular momentum of the disc to j^Rh) as given by 
equation ^ at fixed Vc, then the collapse factor of the disc is 
related to the spin parameter of the halo according to 

aRn^^, (8) 

(Fall & Efstathiou 1980). If the gas conserves its angular 
momentum during collapse, then a disc that forms in a halo 
with a low value of A will collapse by a large factor so ending 
up with a high mean surface mass density; likewise a disc 
that forms in a halo with a high value of A will end up 
with a low mean surface mass density. The large scatter 
in A from tidal torques should therefore be reflected in the 
spread of surface mass densities and surface brightnesses of 
disc systems. The disc luminosity of an exponential disc is 
related to the central surface brightness, la, according to 
Ld = 27r/oa~^, thus from equations |^ and (|), 

/o4''"' oc ^. (9) 

Fig. 1^ shows the distribution of Iqv^'^'^'' for the Byun- 
Freeman sample, compared to the distribution of 1/A^ for 
haloes in the N-body simulations. The mean of the 1/A^ 
distribution has been adjusted to match the mean of the 
observational data points. The observed distribution plot- 
ted in Fig. ^ peaks at about ICQ Lq/pc^ and has a broad 
distribution comparable to, but not quite as broad as the 
distribution of . It is easy to think of reasons why the 
observed distribution might be artificially narrow, for ex- 
ample, the Byun- Freeman sample could be incomplete at 
low and high surface brightness, and perhaps some objects 
with large collapse factors produce spheroidal systems. In- 
deed there is a large literature on observational selection 
effects and the surface brightness distributions of galaxies 
{e.g. Disney 1976; Allen & Shu 1979; Phillipps & Disney 
1986, Dalcanton, Spergel and Summers 1997). Qualitatively, 



however, the above comparisons show that the mean angu- 
lar momenta of disc systems, and their scatter, fit well with 
a picture in which disc galaxies form at recent epochs by 
the collapse of extended gaseous atmospheres within dark 
haloes. Fig. |2|b suggests that disc galaxies form within dark 
haloes with a wide range of A values, i.e. that the angular 
momentum of the halo is not a critical factor in determining 
whether a disc system forms. Within this picture, however, 
the angular momentum of the disc component must be ap- 
proximately conserved during collapse. 

In the following sections, we will investigate this pic- 
ture numerically in some detail, showing that cooling and 
feedback processes are crucial in determining the angular 
momentum of the disc component, as well as in determining 
the fraction of the baryonic gas that can cool. 



3 NUMERICAL TECHNIQUES 

Our cosmological simulations were performed using two 
different self-consistent, 3-dimensional, N-body codes, 
TREESPH (Hernquist & Katz 1989) and GRAPESPH {e.g. 
Steinmetz 1996). These codes are capable of evolving sys- 
tems containing both coUisionless material {i.e. stars and 
dark matter) and gas. 

Smoothed particle hydrodynamics, in which gas is par- 
titioned into fluid elements, is used to model hydrodynam- 
ical interactions between gas particles (Lucy 1977, Gingold 

6 Monaghan 1977). The self-gravitating fluid elements are 
represented as particles that evolve according to Lagrangian 
hydrodynamic conservation laws which take into account lo- 
cal effects arising from pressure gradients and viscosity. The 
equation of state is that for an ideal gas, P = (7— l)p« where 

7 = 5/3, p is the gas density, and u is the speciflc thermal 
energy. As shocks can arise in the gas, an artificial viscos- 
ity is required to maintain numerical stability. The form of 
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artificial viscosity used in both codes is similar to that used 
by Navarro & Steinmetz (1997), i.e. modified by the curl of 
the velocity field to reduce the shear component. 

Radiative cooling and star formation are included in 
the codes. The primordial gas is assumed to have a hydro- 
gen fraction of X = 0.76 and a helium fraction of Y = 0.24 
by weight. CooUng is computed using the ionization, recom- 
bination, and cooling rates from Black (1981, Table 3) with 
the modifications suggested by Cen (1992) for temperatures 
T xi 10^. Radiative cooling docs not proceed below a min- 
imum temperature Tmin = 10, OOOK. The star formation is 
implemented differently in the two codes and will be de- 
scribed in Sections 3.1 and 3.2. 



3.1 TREESPH 

TREESPH computes gravitational forces between parti- 
cles using a hierarchical tree algorithm that groups dis- 
tant particles into nested colls. The potentials of the colls 
are then approximated with a multipole series truncated at 
the quadrupole. When computing the force on a particle, 
TREESPH calculates the ratio of the current cell size to the 
distance between the cell and the particle. For values < 0, 
a specified tolerance parameter, the force from the cell is 
treated as a whole; otherwise, the next cellular subdivision 
is considered. A spline kernel is used to soften the gravita- 
tional potential in order to reduce two-body relaxation. 

In TREESPH, individual smoothing lengths, hi, are cal- 
culated for each gas particle and updated to ensure that 
Nsmooth = 25 — 45 neighbours are contained within a 
sphere of radius 2hi. In addition to individual smoothing 
lengths, individual timesteps are permitted for each parti- 
cle. The largest time step for these TREESPH simulations 
is At = 6.24 X 10® years; the smallest number of timesteps 
in which a particle can reach a = is Ngtep = 2000. Smaller 
timesteps are chosen for particles in dense regions, using the 
criterion aiViAti < etoiEi where Ui and Vi are the accelera- 
tion and velocity of the particle, and etoi = 0.1 determines 
the fractional accuracy of the integrations. Ei is the maxi- 
mum of the mean specific kinetic energy and one sixth of the 
mean specific potential energy of the system (Ewell 1988). 
Particles in the highest density areas are allowed to reduce 
their timesteps by up to a factor of eight. 

Thermal energy is computed from 



ViW{rij,hj)] 



Pi 



(10) 



where the first term is due to adiabatic processes, the sec- 
ond term accounts for the viscosity, and the third term in- 
cludes the remaining nonadiabatic radiative cooling. W is 
the SPH smoothing kernel, which is implemented in the 
gather-scatter formulation of Hernquist & Katz (1989). A 
semi-implicit method of integrating the energy equation is 
employed. An explicit integration using a timestep that is 
half of the shortest SPH particle timestep is performed on 
the adiabatic terms whereas an implicit integration is per- 
formed on the nonadiabatic thermal heating and cooling 
processes. In high density regions, the cooling time step can 
become so short that it is computationally impractical to 
employ it. Thus, in order to avoid numerical errors in the 
integration of the thermal energy, the cooling rate is damped 



so that a particle can lose no more than half its thermal en- 
ergy during a cooling time step (Katz & Gunn 1991). 

Star formation, in which gas particles which meet the 
formation criteria are allowed to form a 'star' wholly and in- 
stantaneously, has been included in TREESPH. The criteria 
for star formation are similar to those of Katz (1992). First, 
gas particles in collapsing regions are found, using V • Vi < 
as the condition for the existence of a convergent flow. Sec- 
ond, the flow is tested for Jeans instability. The particle is 
locally Jeans unstable if the sound crossing time is shorter 
than the dynamical time. Each gas particle in a converging 
flow with 



> 



1 



(11) 



where Ci is the local sound speed, is thus a candidate for 
star formation. However, regions in which the gas density 
has decreased drastically during previous star formation will 
not contract as readily as physically expected due to grav- 
itational softening. Because gas particles in these high star 
formation regions will otherwise be prevented from forming 
stars through inability to meet the Jeans instability crite- 
rion, that criterion is not applied to gas in softened regions 
with 



(12) 



where Ks = 0.89553 is a constant derived from a calculation 
of the Jeans mass (Katz 1992), e, is the softening length, mn 
is the mass of a hydrogen atom, k is the Boltzmann constant, 
and Ti and p; are the gas temperature and density. Third, 
if the cooling time, tc — Ui{dui / dt)~^ , is shorter than the 
local dynamical time, ta = (47rGpf)~'^^^, then gas pressure 
will not inhibit the collapse of the gas particle and it may 
form a star. However, because radiative cooling is cut off 
at Tmin = 10, OOOK, gas particles at low temperatures have 
long cooling times. In order to avoid an artificial cessation of 
star formation in low temperature regions, gas particles with 
T < 30, OOOK which otherwise meet the criteria are allowed 
to form stars with a cooling time set to the dynamical time 
tj, =td. The local rate of star formation is 

dp* ~dpg 
It ~ 



dt 



tc 



(13) 



where c, = 0.1 is a dimensionless star formation parameter. 
The probability that a gas particle will form a star in time 
At is 



(14) 



A random number is generated for each star candidate; if it 
is less than p, , the gas particle loses its SPH characteristics 
and becomes a coUisionless star particle. 

For these cosmological simulations, TREESPH evolves 
the haloes from the initial redshift to « = using physical 
coordinates with isolated boundary conditions. 

3.2 GRAPESPH 

The version of GRAPESPH we employ has evolved from 
the TREESPH code of Navarro & White (1993; hereafter 
NW). Rather than using the mutual nearest-neighbour bi- 
nary tree (Benz et al. 1990) of the original implementation 
in order to calculate gravitational forces and find neighbours 
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for the SPH force evaluation, the special-purpose GRAPE- 
3Af hardware (Sugimoto et al. 1990) has been harnessed to 
perform these tasks. This machine carries out a parallelised 
direct summation over all particles and very rapidly returns 
the forces and neighbour lists. All hydrodynamical forces 
are then calculated on the host workstation. The adapta- 
tions which convert a TREESPH code into a GRAPESPH 
code have been described in detail by Steinmetz (1996). 

Both individual timesteps and individual gas particle 
smoothing lengths were used. As the GRAPE only returns 
'gather' neighbours, that is those within distance hi of gas 
particle i, the neighbour search radius is taken to be \/l-5 
times the smoothing length. This should ensure that the 
actual number of mutual neighbours, which satisfy hi + hj < 
twice the particle separation, lies within the range 25 — 45. 

The gravitational forces are calculated assuming a 
Plummer softening law. For each particle, the total accel- 
eration due to gravity is found in two steps. First, the ac- 
celeration arising from all gas particles is found, and then 
that resulting from dark matter and stars is incremented. A 
different gravitational softening parameter is used in these 
two separate force evaluations. 

A simplified version of the star formation algorithm de- 
scribed by NW was adopted. Any gas particle in a suffi- 
ciently high density, collapsing region is considered to be a 
potential star-forming particle. These two constraints can 
be expressed as 



Pi > Pc: 

and 



= 7 X 10" 



kg m 



V.V; < 0. 

The local dynamical time is then defined by 



'^dyn — 



16Gpi ■ 



(15) 



(16) 



(17) 



If conditions (^5|) and (^^ remain satisfied for the ensuing 
dynamical time then the gas particle converts entirely to a 
star particle. Should either of these constraints be violated 
during this probationary period then the gas particle is no 
longer considered to be star-forming and it remains gaseous. 

3.3 Contrasting the Codes 

The softening of the gravitational interaction between two 
particles differs in the two codes. In GRAPESPEI, the po- 
tential is of the form corresponding to a Plummer density 
profile, $ oc {r^ + e^)~^^^. However, the acceleration cal- 
culated using a Plummer model softening converges to the 
Kepler value relatively slowly. In codes like TREESPH, in 
which distant particles are assigned to cells and their po- 
tentials are represented by quadrupole expansions, the cells 
are assumed to be point-like rather than extended. Thus, a 
softening in which the acceleration quickly converges to the 
Kepler form is advantageous. In TREESPH, the gravita- 
tional potential is softened using the spherically symmetric 
spline kernel suggested by Monaghan & Lattanzio (1985) be- 
cause it has compact support and is identical to the Kepler 
form for r > 2e. For the effective softening to be similar in 
both codes, the TREESPH softening lengths must have val- 
ues approximately twice those of the GRAPESPH softening 
lengths. 



The star formation algorithms employed within the two 
codes derive from Katz (1992) and NW. They are simi- 
lar in requiring that a gas particle be able to form a star 
only if it is in a converging fiow. However, other criteria 
differ. This version of GRAPESPH, unlike TREESPH and 
NW, does not use the Jeans instability criterion for creating 
star particles nor does it define a star formation rate. Like 
NW, GRAPESPH defines a critical density, above which 
the cooling time is less than the dynamical time. The two 
GRAPESPH criteria are applied to each particle when it 
is advanced. The criteria are monitored over the dynamical 
time for each 'successful' gas particle after which, if they are 
not violated, a star is formed. In general, the delay is ^ 10^ 
years. In TREESPH, the criteria are checked at each sys- 
tem timestep of At = 6.24 x 10® years; a random number of 
gas particles which meet the star-forming criteria are chosen 
to form stars immediately. The effect of these differences is 
that the GRAPESPH simulations produce more stars be- 
tween z = 1 and 0.8 while, afterwards, the star formation 
rate in the TREESPH simulations is higher. By z = 0, the 
number of GRAPESPH stars is some 75—90% of the number 
of TREESPH stars. This difference is due to the additional 
TREESPH formation criteria which allow stars to continue 
to form in regions of high star formation where the gas den- 
sity has decreased due to the transformation of gas particles 
to stellar particles. Qualitatively, the GRAPESPH criteria 
act to mimic larger supernova feedback effects in regions of 
high star formation. 



4 INITIAL CONDITIONS 

The initial conditions for our simulations are generated us- 
ing a nested hierarchy algorithm similar to that described 
by Katz & White (1993). The idea is to select haloes from 
the dissipationless CDM simulation described in Section ^ 
and then to re-simulate them at higher mass resolution us- 
ing the SPH codes of the previous section. This technique 
has been used by Navarro, Frenk & White (1995b) and Eke, 
Navarro & Frenk (1998), in SPH simulations of clusters of 
galaxies and by Navarro & Benz (1991), Quinn, Katz & Ef- 
stathiou (1996), Navarro & Steinmetz (1997) and others in 
SPH simulations of galaxy formation. 



4.1 Selection of the Haloes 

Haloes with a density contrast A ~ 200 were located at 
z = in the N-body simulation using the group finding al- 
gorithm described in Section |^. The haloes were ordered by 
circular speed. We then visually examined plots of the halo 
particles at various output times and selected five examples 
for simulation with the SPH codes. The haloes were chosen 
to have a range of circular speeds between 300 kms~^ and 
150 kms~^, to be far from much more massive haloes, and 
not to have merged with a comparable mass system since 
z = 1. Parameters for the 5 chosen haloes are given in Table 
1. The first condition selects objects with circular speeds 
in the range observed for normal luminous disc systems. 
The second criterion eliminates satellite systems near the 
outer parts of massive systems, for which our tree codes and 
hierarchical initial conditions algorithm are ill-suited. The 
third criterion eliminates haloes that undergo major merger 
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Figure 3. The evolution of two haloes from the dissipationless N-body simulation described in Section ^ x-y projections are shown of 
the particles located within a box of comoving length 4 Mpc centred on the dominant subclump at 2 = 2, 2 = 1 and 2 = 0. The upper 
panel shows the evolution of Halo 2 listed in Table 1, with a final circular speed of Vc = 249 kms~^. The lower panel shows the evolution 
of a halo with a nearly identical final circular speed, Vc = 247 kms~^, that we did not simulate with the SPH codes because it evolves 
by the merger of two comparable mass components at 2 < 1. 



events at low redshifts and so preferentially selects regular 
objects that grow by the accretion of smaUer systems at low 
redshifts. Nevertheless, even with this selection criterion, the 
typical halo mass within an overdensity of A = 200 nearly 
doubles between a redshift of 1 and (see Table 1). 

Table 1: Parameters for the five haloes 
chosen for SPH simulation 



the box used to gener ate h igh resolution initial conditions 



Halo 



(km/s) 



Mh X 1O'M^0 

2 = 2 = 1 



\h Lhr (Mpc) 



1 


265 


70.5 


35.2 


0.050 


6.00 


2 


249 


57.7 


38.8 


0.020 


6.08 


3 


176 


25.1 


14.2 


0.031 


4.92 


4 


160 


21.5 


11.8 


0.040 


5.64 


5 


157 


20.1 


16.7 


0.032 


4.84 



Notes to Table 1: The first column gives the halo number, 
used throughout this paper. The second column gives the 
circular speed computed from equation ^. The third column 
gives the mass of the halo at 2 = 0, and the fourth column 
gives the mass of the dominant subclump at 2 = 1 identified 
by running the group finding algorithm with a link param- 
eter b — 0.2. The fifth column gives the A parameter of the 
halo at 2 = 0. The final column gives the comoving size of 



as described in Section 4.2 



The upper panel in Fig. |^ shows an x-y projection of 
Halo 2 at 2 = 2, 2 = 1 and 2 = 0. The lower panel shows 
the evolution of a halo with almost identical circular speed 
that we decided not to simulate because it evolves by the 
merger of two comparable mass subunits at 2 < 1. The 
selection criteria that we have imposed will bias statistics 
of the frequency of disc galaxy formation in CDM models, 
and this must be borne in mind when interpreting some of 
the results presented in the next two Sections. In particular, 
the last column of Table 1 shows that our selection criteria 
bias against haloes with high values of A. This is because 
haloes with high A often contain two distinct subclumps, or 
are near much more massive systems. 



4.2 Generation of High Resolution Initial 
Conditions 

The particles within a sphere of radius 400 kpc around the 
centre of mass of the halo at 2 = were traced back to 
the initial conditions. The smallest cubical box containing 
all of these particles was located, defining the region to be 
resampled at higher mass resolution. The comoving box sizes 
of the high resolution regions for each halo are listed in the 
final column of Table 1; typically a box of comoving size 
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5-6 Mpc is required. Having located the size and location 
of the high resolution region, initial conditions for the SPH 
simulations were generated as follows. 
[1] Distribute N^^ particles on a regular lattice within the 
high resolution box of size Lhr- 

[2] Generate displacements of these particle positions using 
the algorithm of Efstathiou et al. (1985) with the identical 
amplitudes and phases used to generate the initial conditions 
of the low resolution simulation, but with the initial power 
spectrum sharply truncated at the Nyquist frequency of the 
low resolution particle grid kn = nNp/Lbax- 
[3] Generate additional displacements to add in the missing 
small-scale power between the Nyquist frequencies of the 
low resolution particle grid and the high resolution parti- 
cle grid kn' = -KNhr/Lhr from the appropriate CDM power 
spectrum. 

[4] Add M nested layers of massive particles to represent the 
tidal field of the surrounding matter. The particles in these 
outer layers are arranged on lattices such that the lattice 
spacing increases by factors of two outwards from the high 
resolution region and particle masses increase by factors of 
eight, as described by Katz &: White (1993, see their Fig. 
1). Five layers were chosen to match exactly the total mass, 
Afbox, contained in the low resolution N-body simulation. 
[5] Generate displacements for these outer particles as in [2] 
above, but with the power spectrum sharply truncated below 
the Nyquist frequency of each layer. This last step prevents 
aliasing of small-scale power from affecting the representa- 
tion of the tidal field. 

For the simulations described in the next section, we 
used a high resolution region represented by 27^ particles. 
The outermost of the five low resolution layers of dark mat- 
ter particles was stripped off each halo, because tests re- 
vealed that it had a negligible effect on the evolution of the 
particles within the high resolution region. The total mass, 
A/tot ! of the four remaining layers and the interior high reso- 
lution particles was calculated and the new outermost layer 
was truncated by removing all particles within radius 



Mu 



1/3 



(18) 



2 VAfbo 

from the centre of the simulation. This produces a spheri- 
cal distribution of the outermost high mass particles. Each 
dark matter particle in the high resolution layer was re- 
duced in mass by 10% and overlaid with a gas particle of 
mass rugas = Q.lmdark to produce a high resolution region 
containing gas. The masses of dark matter and gas parti- 
cles in the simulations range from 3.6 — 7.1 x 10* Mq and 
4.0 - 7.9 X 10'^ M0, respectively. 

The total number of particles is Ntot = 41136, with 
the low resolution layers having A'^i = 1178, N2 ~ 386, 
N3 = 152, and N4 = 54. In the TREESPH simulations, the 
gas particles were offset in position from the dark matter 
particles by a fraction of a softening length to avoid numeri- 
cal difficulties in constructing a tree containing all particles. 



5 NUMERICAL SIMULATIONS 

5.1 Parameters of The Simulations 

Each of the five haloes of Table 1 was simulated at high 
mass resolution with TREESPH and GRAPESPH. As we 



will show in subsequent sections, the softening of the grav- 
itational force law has a large effect on the calculations. In 
the TREESPH simulations, the high resolution dark mat- 
ter particles are assigned gravitational softening lengths of 
Edarfe = 10-4 kpc. The gas and star particles in TREESPH 
have egas = estars = 2.6 kpc to allow better resolution of the 
overdense regions in which the galaxies form. Larger soften- 
ing lengths are adopted for the low resolution layers of high 

mass particles according to eiayer = {miayer/rndark)^^'''^dark- 

In the GRAPESPH siniulatious. the gravitational forces 
are computed assuming a Plummer force-law with soften- 
ing lengths of €dark,star = 5 kpc for the dark matter and 
star particles and egas = 1 kpc. With these choices, the soft- 
ened force laws for the gas and high resolution dark matter 
particles are similar in the two codes. However, the softening 
adopted for the stars in GRAPESPH is about twice as large 
Eis the softening adopted for the star particles in TREESPH. 
This is because in the GRAPESPH code, gravitational forces 
for all coUisionless particles (dark matter and stars) are com- 
puted in one single parallelized direct summation and thus 
have the same softening. 

Both codes employ variable timesteps. Particles in the 
TREESPH runs completed the evolution from the initial 
redshift of « = 7.4 to a = in 2000 to 16000 timesteps and 
a typical simulation with 41136 particles required approx- 
imately 100 hours to run on an DEC Alpha 4100 5/300. 
Evolution of the particles in the GRAPESPH simulations 
required typically ~ 100 - 40000 steps and ^ 45 - 65 CPU 
hours to run to completion from z = 7.4 to 2; = on the 
GRAPE machine at Edinburgh. 

In most previous simulations of this type, the gas has 
been allowed to cool radiatively with no other energy input 
except that provided by gravitational shock-heating and, in 
some simulations, a uniform photoionizing background radi- 
ation (see e.g. Quinn et al. 1996; Navarro & Steinmetz 1997). 
In some simulations, attempts have been made to mimic en- 
ergy injection into the intergalactic medium by supernovae 
(Katz 1992; NW). Broadly speaking, these have involved ei- 
ther injecting energy directly as heat into the gas, or into 
bulk motion so reversing the flow of gas into high density re- 
gions. As discussed by NW, neither approach is satisfactory. 
Energy supplied as heat into high density gas is radiated 
away so efficiently that it has very little feedback effect. On 
the other hand, it is easy to reverse the flow of gas by in- 
jecting energy into bulk motions, but the feedback effects 
arc then extremely sensitive to the specific algorithm and 
parameters that are used. We believe that a more realis- 
tic description of feedback requires the ability to model a 
multiphase gas component, in which cool high density gas 
can coexist spatially with an outflowing hot lower density 
gas component. Modeling such a nnilti-phase gas medium 
would require fundamental modifications to the SPH codes 
described here, and so is well beyond the scope of this pa- 
per. Some provisional attempts along these lines, using a 
Eulerian gas dynamic code, are described by Yepes et al. 
(1997). 

Nevertheless, as described in the Introduction, our the- 
sis in this paper is that feedback processes are critical in 
determining the angular momenta of disc systems. The key 
process that we wish to model is the late infall of gas into a 
dark matter halo. The numerical simulations of Navaxro & 
Steinmetz (1997) described in the introduction show that, in 
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the absence of any feedback, gas collapses at high redshifts 
into dark matter subclumps that later merge, losing most of 
their orbital and internal angular momentum to outer halo 
particles. If, however, feedback processes can prevent some 
significant fraction of the gas from collapsing at high red- 
shifts, then there is a possibility that this residual gas can 
collapse within a more uniform and slowly evolving halo to 
form a disc system whilst conserving its angular momentum. 
Such a formation scheme for disc systems seems essential in 
hierarchical clustering theories, since the discussion of Sec- 
tion ^ shows that the angular momenta of disc systems can 
only be understood if the gas approximately conserves its 
angular momentum during collapse. 

In this paper, we have chosen a particularly simple 
scheme to model feedback processes. We evolve the gas com- 
ponent adiabatically until a specified redshift Zcooi, at which 
point we switch on radiative cooling with a cooling rate A 
as described in Section ^. For most of the runs, we adopted 
Zcooi = 1. Some runs were done with Zcooi = 4 and Zcooi = 0.6 
to illustrate the sensitivity of angular momentum evolution 
to the epoch of gas collapse. Two simulations were done in 
which the cooling rate A was suppressed by a factor, 

g{z) = {l + zrei^p{~{z/z,f), a = -1.5, (19) 

where Zc is a parameter. The parameter a was fixed to 
a = —1.5 by an analytical calculation of cooling within an 
isothermal halo so that discs would grow roughly linearly 
with time. Simulations were done with Zc = 1, almost com- 
pletely eliminating high redshift cooling, and with Zc = 4. 
Table 2 gives a complete list of the SPH simulations that 
were run for this paper. 

5.2 Comparison of Codes 

The simulations described here are at the limit of what is 
possible with present-day computers, thus it is important to 
understand the limitations imposed by the computational 
scheme, such as limited mass resolution, gravitational soft- 
ening etc. In this subsection, we thus describe in some detail 
the differences in the evolution of two haloes (2 and 3) run 
with the two codes. 

5.2.1 Morphologies 

Fig. ^ shows projections of the dark matter and gas distri- 
butions ai z = 1 evolved with cooling suppressed (Runs 4, 
5, 11, and 12). The left-hand panels show the particle distri- 
butions evolved with TREESPH and the right-hand panels 
those evolved with GRAPESPH. Although the two codes 
produce dark matter and gas distributions that are closely 
similar, we were surprised to see some differences. Perhaps 
the most striking variations can be seen in the relative po- 
sitions and morphologies of the four dark matter concentra- 
tions to the left of Halo 3. The structure of the central halo 
itself also appears to be slightly different in the two runs. 

To explore these morphological differences, and deter- 
mine whether they were caused by differences in the treat- 
ment of the gas component, we re-evolved Halo 3 with the 
two codes from z = 7.4 to z = 1 using only dark matter 
particles. First, the haloes were run with the same compu- 
tational parameters as Runs 11 and 12. Projections at z = 1 



Table 2: List of Simulations 



Run Halo Code ^cooi Zc 

1 
2 
3 
4 
5 
6 
7 
8 
9 

10 
11 
12 
13 
14 
15 
16 
17 
18 



Notes to Table 2: The first column lists the run number, the 
second gives the halo number as in Table 1 and the third 
column lists code used to run the simulation. The fourth 
column lists the redshift Zcooi at which radiative cooling was 
switched on. Where this is blank, we suppressed the cooling 
rates by the factor given in equation [13 with parameter Zc 
listed in the fifth column. 

are shown in panels (a) and (b) of Fig. ^, which can be com- 
pared directly with the dark matter plots shown in Fig. ^ 
The differences between the particle positions computed 
from the two codes are almost identical to those seen in 
Fig. 1^, implying that they do not arise merely from differ- 
ences in the treatment of gas. Rather, the differences must 
arise from errors in the gravitational forces or the time- 
integration of the equations of motion. We therefore reran 
these dark matter only runs using different code parame- 
ters. For TREESPH, we improved the accuracy of the forces 
by changing the force-tolerance parameter from the default 
value 6* = 0.8 to = 0.5. In GRAPESPH, the force calcula- 
tion is hardwired into the GRAPE with, typically, a ~ 2% 
dispersion in the forces between particle pairs (see Okumura 
et al. 1993). Since the force accuracy cannot be adjusted in 
the GRAPESPH code, we reran the simulation employing 
three times as many timesteps. Snapshots from these two 
runs are shown in panels (c) and (d) of Fig. |^, and show 
clearly that force errors in TREESPH are the main cause of 
the discrepancies. For many purposes these small errors are 
unimportant, but they can lead to different merger histories 
and to differences in the radial profiles of some quantities 
sensitive to the presence of substructure, e.g. plots of the 
specific angular momentum with radius (see Figures 8 and 
12 below). 

The stellar condensations that form in Runs 4 and 5 
(Halo 2) and Runs 11 and 12 (Halo 3) are shown in Fig. ^ 
These 'galaxies' have been rotated so that each panel shows 
a projection along two of the principal axes; the TREESPH 
X — Y and X — Z projections are shown on the left, and the 
corresponding GRAPESPH X — Y and X — Z projections 
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Figure 5. Four simulations of Halo 3, using dark matter only, 
plotted at 2 = 1. Panels (a) and (b) show TREESPH and 
GRAPESPH using the default computational parameters. These 
plots are almost identical to the dark matter distributions plotted 
in Fig. ^ and show distinct differences between the two codes, for 
example, in the location of the four clumps to the left of the main 
central mass concentration. Panels (c) and (d) show a TREESPH 
simulation employing a more accurate force computation, and a 
GRAPESPH simulation using smaller timesteps. Panels (c) and 
(d) are almost identical, and similar to the default GRAPESPH 
simulation, showing that force errors in TREESPH cause most of 
the discrepancies. 



Figure 4. Projections of the dark matter and gas components 
at z = 1 for two haloes (Runs 4, 5, 11, and 12) evolved with the 
TREESPH and GRAPESPH codes. 



are shown on the right. Halo 2 evidently forms a disc sys- 
tem in both simulations. Similar numbers of stars are formed 
with the two codes, but there are noticeable differences at- 
tributable mainly to the much larger stellar softening length 
used in GRAPESPH. The TREESPH galaxy is clearly more 
concentrated, and has a much thinner disc. 

Much larger differences are evident in the evolution of 
the Halo 3 stellar systems, shown in the lower panels of 
Fig. ^. The TREESPH galaxy shows disc-like structure at 
both redshifts. At 2 « 0.5, the main component in the 
GRAPESPH simulation is also disc-like, but there is a satel- 
lite system containing about 25% of the mass of the main 
galaxy which is not apparent in the TREESPH panel. The 
satellite system is present in the TREESPH simulation, but 



the slight difference in output redshifts (0.49 as opposed to 
0.46 for GRAPESPH) coupled with the high relative speed 
of the intruder conspire to move it outside the region plot- 
ted in the figure. In both runs the incoming satellite makes 
a first fly-by of the main disc at z ~ 0.4, with the even- 
tual merger occurring at z ~ 0.15. However, because of the 
more extended structure of the GRAPESPH central stellar 
disc that results from the larger softening, the interaction 
and subsequent accretion of the satellite system is more dis- 
ruptive in this case. As a consequence the merger event has 
a preferentially destructive effect on the GRAPESPH disc, 
leaving a spheroidal stellar object at 2; = 0, whereas the 
TREESPH disc survives largely intact. Referring back to 
Fig. ^ it is the subclump directly below the central clump 
that is responsible for this collision (the systems to the left 
do not fall into the centre before z = 0). This example shows 
the sensitivity of the final remnant morphologies to the stel- 
lar softening. 
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Figure 6. Stars in the Halo 2 (Runs 4 and 5) and Halo 3 (Runs 11 and 12) simulations at z f» 0.5 and 2 = 0. Lefthand panels sliow the 
TREESPH simulations projected along the principal axes X ~Y and X — Z and the righthand panels show the GRAPESPH simulations 
in the same projections. 



5.2.2 Star Formation 

We define a 'virial' radius, Rvir{z), for each halo by placing 
a sphere at the centre of mass within which the overdensity 
is equal to A = 178. The number of star particles in our 
simulations within the virial radius at z — lies within 
the range Nstar = 1700 - 4000. In Fig. the stellar mass 
within a virial radius for the Halo 2 and 3 galaxies is shown 
over the redshift range z — 1 to z — 0. The two codes 
produce a comparable number of stars before z = 0.8, after 
which the TREESPIf haloes accumulates stellar mass more 
rapidly. In both codes, stars are accumulated at a faster 
rate in Halo 2 than in Halo 3. The Halo 2 galaxy has a 
final stellar mass of 2.4 x IO^Mq in GRAPESPH and 2.7 x 
10" M© in TREESPH. For Halo 3, the final stellar mass is 



1.17 X 1O"M0 in the GRAPESPH run and 1.45 x IO^Mq 
in the TREESPH run. Halo 2 forms the more massive stellar 
disc, by about a factor of two. 

It is interesting to compare these numbers with the stel- 
lar masses of L* galaxies. The luminosity function of Efs- 
tathiou, Ellis and Peterson (1988) gives L% = 5.1 x lO^^L© 
for h = 0.5. Assuming that the mass-to-light ratio of typical 
discs is {M/L)b ~ 3h in solar units, the stellar mass of a 
typical L* disc system is about M*disc ~ 8 x 10^°Mq, i.e. 
somewhat smaller than the masses of the stellar discs that 
form in Halos 2 and 3. The agreement is close enough, how- 
ever, that it suggests that the basic physical picture of late 
infalling gas may account for the masses of disc systems and 
the origin of the TuUy-Fisher relation (equation W. 

Fig. ^ also shows the stellar mass within a virial radius 
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Figure 7. Stellar mass within a virial radius for the Halo 2 and 
3 galaxies between redshifts z = I and 2 = 0. Dotted lines are for 
TREESPH simulations of Halo 2 using the cooling suppression 
formula, g{z) with Zc = 1 (bottom line) and with Zc = 4 (top 
line). 
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for the two TREESPH simulations of Halo 2 using the cool- 
ing suppression formula, g(z) of equation (|l^). The bottom 
dotted line is for a cut-off redshift of Zc — 1 and the top dot- 
ted line is for Zc = 4. In both these simulations, radiative 
cooling and star formation are turned on at 2; = 4. However, 
star formation is clearly nearly suppressed in the Zc — 1 
model until z — 1, after which stars continue to form slowly 
until z — 0. In the Zc — 4 model, nearly half the final total 
mass of stars has formed by z = 1, but by z = the galaxy 
consists of slightly fewer stars than in the simulation with 
cooling suppressed until Zcooi ~ 1- 



5.2.3 Angular Momentum Profiles 

The cumulative specific angular momentum profiles of the 
stellar components in Haloes 2 and 3 are shown in Fig. ^ 
The total J/AI within a sphere is calculated for 90 equally- 
spaced radii. The specific angular momenta measured in the 
TREESPH simulations are about 50% higher than in the 
GRAPESPH simulations. A small companion at a distance 
r « 170 kpc appears as a stepwise increase in specific angular 
momentum of Halo 3 in the TREESPH run; this difference 
to the GRAPESPH run is a consequence of the small differ- 



ences in the dark matter evolution discussed in Section f>.2 
In Section 



we will show that angular momentum evo- 
lution of the stellar and gaseous components can vary by 
an order of magnitude or more, depending on what assump- 
tions are made about the physical parameters of the gas 
{e.g. suppression of cooling). These differences in J/M aris- 
ing from the input physics are therefore much greater than 
the differences between the codes seen in Fig. H. 

The most obvious differences in Fig. ^ are in the spe- 
cific angular momenta in the central regions, r ^ 20 kpc. On 



Figure 8. Cumulative specific angular momentum profile of the 
stellar component out to 200 kpc for the Halo 2 and Halo 3 galax- 
ies at 2 = 0. 



these small scales, the J/M profile for the GRAPESPH runs 
declines much more steeply than in the TREESPH runs. Ev- 
idently, the angular momentum profiles decline rapidly on 
scales smaller than a few softening lengths, and this effect 
is much more pronounced in GRAPESPH because of the 
muc h larg er stellar softening length. As described in Sec- 



tion 5.2.1 



the large softening length of GRAPESPH leads 
to more diffuse stellar systems, and clearly one cannot ex- 
pect these simulations to describe accurately the mass and 
velocity profiles of real disc systems. However, global prop- 
erties, such as the total stellar mass and angular momentum 
are much less sensitive to softening. 



5.2.4 Rotation velocities and velocity dispersions 

Fig. P shows the projected velocity profiles out to 19 kpc 
of the Halo 2 and 3 galaxies. The projected velocities are 
calculated by distributing the particles onto a Cartesian grid 
in which X, Y, and Z are chosen to be along the major, 
intermediate, and minor principal axes, respectively. Each 
grid cell has dimensions of A/ = 2 kpc. A slit with a width 
of two cells was laid parallel to the axis; results are averaged 
over the two cells. In Fig. ^, the mean rotation speed of 
the stars (solid lines) and the projected velocity dispersion 
(dashed lines) are plotted against major axis distance X, 
i.e. as would be measured for an edge-on galaxy. 

Clear evidence for systematic rotation is seen in both 
TREESPH runs. In the GRAPESPH runs, the Halo 2 galaxy 
is rotating with a peak amplitude of about 80 kms~^. The 
difference with TREESPH is, again, primarily a consequence 
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Figure 9. Projected major-axis velocity profiles out to 19 kpc 
for tlie Halo 2 and 3 stellar systems at 2 = 0. Solid lines show 
the rotation velocity, Vr \ dashed lines show the projected velocity 
dispersion, u. 



of the larger GRAPESPH softening length, which inhibits 
the formation of a compact, rapidly rotating stellar system. 
Little rotation is seen in the GRAPESPH Halo 3 stars, as 
might be expected from the lack of any disc in Fig. ^. In 
comparison, the TREESPH Halo 3 disc rotates with a max- 
imum rotation speed of Vr ~ 100 km/s at X « 3 kpc. In 
both codes, the projected velocity dispersions he in the range 
a ~ 100 — 200 km/s, and are roughly constant with radius. 
The velocity dispersions are comparable to the mean rota- 
tion speed in the TREESPH code; thus, neither code suc- 
ceeds in producing a cool, centrifugally supported, stellar 
disc. This is almost certainly caused by two-body heating 
(see e.g. Steinmetz & White 1997) which, for discs com- 
posed of a few thousand particles, can heat a cold disc to to 
Vr/<y ~ 1 within a few rotation periods. This is consistent 
with the time-evolution of a which increases by a factor of 
« 2 between redshift z = 1 and 0. 



5.2.5 Stellar Mass Profiles 

Fig. |l^ shows the cumulative stellar mass profile for the Halo 
2 and 3 galaxies out to a radius of 25 kpc. This gives an idea 
of the concentration of the final stellar remnants and of the 
large differences caused by the softening in GRAPESPH. 
Even in the TREESPH code, the half-mass radii of the stel- 
lar systems are ~ 3 kpc, about equal to stellar softening. 
This shows that it is possible to form rapidly rotating discs 
with similar scale lengths to real disc galaxies, though higher 
resolution calculations are required to model their mass and 
velocity profiles. 
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Figure 10. Cumulative stellar mass profile of the Halo 2 and 3 
galaxies out to a radius of 25 kpc. 



5.2.6 Summary 

The simulation of disc galaxy formation, including star for- 
mation, is a complex problem and it is often difficult to as- 
sess the accuracy of any particular computational technique. 
The purpose of this subsection has been to compare the re- 
sults of evolving identical conditions with two independent 
SPH codes with the aim of finding which properties, if any, 
are sensitive to numerical methods. To summarise the find- 
ings, the dark matter and gas evolution is broadly similar. 
However, in one of the 5 haloes (Halo 3), the small differ- 
ences that do exist produce a significantly different stellar 
remnant. The overall morphologies of the stellar systems 
that form in the other four haloes are more similar. Al- 
though the star formation algorithms are different in the 
two codes, they give similar star formation rates and final 
stellar masses. This is not too surprising because the algo- 
rithms are designed to convert high density cool gas into 
stars on a timescale of order the dynamical time, which is 
much shorter than the Hubble time when most of the stars 
form in these simulations. The star formation rate is there- 
fore governed mainly by the rate at which gas infalls towards 
the halo centres, and this is modeled in a similar way in both 
codes. 

The main differences between the codes are attributable 
to the larger stellar softening employed in the GRAPESPH 
simulations. This results in stellar systems that are more 
extended and have lower rotation velocities than those that 
form in the TREESPH simulations. Even with TREESPH, 
the stellar softening of 2.6 kpc is barely adequate to resolve 
the internal structure of the discs that form. Despite these 
differences, the global properties of the stellar remnants, e.g. 
their masses and angular momenta, are in reasonably good 
agreement. In the next Section we will investigate, using 
both codes, how these global properties depend on the input 
physics. 
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Figure 11. Evolution of the mean specific angular momentum 
within a virial radius of the dark matter (top) and gas (bottom) of 
the five simulations run with ^cool = 1- Here and in all succeeding 
figures, J/M is in kpc km s~^. 



5.3 Angular Momentum Evolution 

We first discuss the evolution of angular momentum in the 
simulations in which cooling is suppressed until Zcooi = 1- 
We then discuss how varying Zcooi alters the conclusions. 

The accumulation of angular momentum by dark mat- 
ter haloes and their associated gas proceeds smoothly dur- 
ing the expansion phase before turnaround, growing roughly 
as t^^^ (White 1984, Barnes and Efstathiou 1987). After 
collapse, halo angular momentum growth slows and then 
evolves by abrupt jumps caused by discrete merger events. 
Since the haloes are all of a similar size, their total specific 
angular momentum within the virial radius, Rvir{z), was 
averaged in order to give a statistical representation of the 
growth of J/M with redshift. The evolutions of the mean 
specific angular momentum for the gas and dark matter are 
shown in Fig. 0. The TREESPH and GRAPESPH resuhs 
are plotted as solid and dashed lines respectively, and are 
almost identical. The specific angular momenta of the dark 
matter and gas have very similar values, but at late times 
the specific angular momentum of the gas component ex- 
ceeds that of the dark matter. This is because J/AI of both 
the components increases with radius. Consequently, when 
the central, low J/M, gas turns into stars, the specific angu- 
lar momentum of the remaining gas within the virial radius 
is enhanced. 



Fig. |l2| shows the radial specific angular momentum 
profiles of the stellar components in each halo for both 
TREESPH and GRAPESPH simulations at three different 
epochs. The redshifts for TREESPH results are z = 0.63 
(dotted lines), 0.30 (dashed lines), and 0.0 (sohd lines); for 
GRAPESPH, the first redshift is z = 0.68. It can be inferred 
from the large number of almost horizontal lines with the 
occasional vertical jump, that the dominant stellar concen- 
trations are usually well isolated from other star particles. 
The two most disc-like final stellar components, belonging 
to Haloes 2 and 5, evolve relatively steadily from z — 0.6 
to the present, whilst the remaining haloes suffer more in- 
teractions with other stellar clumps. The differences in the 
positions and sizes of the vertical jumps in J/A'I in identical 
halos run with TREESPH and GRAPESPH are caused by 
the evolutionary differences identified in Section 5.2.1 that 
lead to sub-condensations in each code with different orbits. 

The morphology and specific arrgular momenta of the 
final stellar objects appear to be sensitive to the merger 
histories undergone by the parent haloes. This point has 
been investigated quantitatively by resimulating two haloes 
with the cooling switched on at Zcooi = 4 and Zcooi = 0.6, in 
addition to the usual Zcooi = 1- The evolution of the stellar 
J/M measured within the central 20 kpc is shown in Fig. |l3[ 
The top frame shows GRAPESPH resuhs for Halo 2; the 
bottom frame shows TREESPH results for Halo 5. For Halo 
2, it is clear that a significant increase in the final specific 
angular momentum occurs when the gas is prevented from 
cooling at early times. The same result is reproduced by the 
Halo 5 simulations, but only for the change of Zcooi from 4 
to 1. This halo evolves very little after z = 1 (see Table 1); 
the result of suppressing cooling until Zcooi = 0.6 is to form 
a very similar stellar disc as in the ^cooi ~ 1 simulation, but 
at a slightly later time. From 2: = 1 to the present. Halo 5 
increases its virial mass by a factor of ~ 1.5 whereas Halo 
2 grows by a factor of ~ 1.2. The other three haloes nearly 
double in mass during this period. 

Fig. ^ shows the specific angular momenta of all the 
three components (stars at 20 kpc, gas and dark matter 
at the virial radius) for both TREESPH (filled symbols) 
and GRAPESPH (open symbols) haloes at redshift zero. 
The circular speed of the haloes is computed from equa- 
tion (1) with Rc set to the virial radius R^ir- Concentrating 
on the stellar components (plotted as 4— pointed stars for 
Zcooi ~ 4, 5— pointed stars for ^cooi = 1 and 6— pointed stars 
for Zcooi = 0.6), it can be seen that, even allowing for the 
differences between the two codes, there is a tendency for 
the more slowly evolving haloes (numbers 2 and 5) to give 
rise to more rapidly rotating stellar objects, provided cooling 
is suppressed at high redshift. If cooling is allowed to begin at 
Zcooi ~ 4, the stellar systems that form in these haloes expe- 
rience the 'angular-momentum catastrophe' seen in previous 
numerical simulations. The crosses in Fig. |l^ show the spe- 
cific angular momenta of the stars for the Halo 2 TREESPH 
simulations in which cooling was suppressed using the em- 
pirical formula g{z) of equation (p^). The lower cross is for 
Zc = 4 and the upper cross is for Zc — 1. The specific an- 
gular momentum of the Zc — 1 simulation is an order of 
magnitude larger than the Zc — 4 run, though it is not as 
high as in the simulation with cooling abruptly truncated 
at Zcooi ~ 1- Nevertheless, these simulations demonstrate 
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Figure 12. Radial specific angular momentum profiles of the stellar components in eacii halo for the simulations with Zcool = 1- Redshifts 
for TREESPH results are z = 0.63 (dotted lines), 0.30 (dashed lines), and 0.0 (solid lines); redshifts for GRAPESPH results are similar 
except that the first redshift is 2 = 0.68. 
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Figure 13. Evolution of the stellar J/M measured within the 
central 20kpc for GRAPESPH Halo 2 and TREESPH Halo 5. 



again the sensitivity of angular momentum evolution to the 
cooling history. 

The solid line in Fig. |l^ shows equation ^ The specific 
angular momenta of the dark matter haloes all lie below 
this line. This is expected because, as described in Section 
4.1, our selection criteria bias against haloes with final spin- 



parameters larger than A — 0.05. All of our haloes therefore 
have spin parameters less than or equal to the median value 
of A measured in N-body simulations (see Figure 1) and 
hence lie below the line plotted in Fig. [141. 

The gas components generally have higher final spe- 
cific angular momenta than their dark matter haloes. As 
explained at the start of this subsection, this is because the 
lower-angular momentum gas in the central regions of the 
haloes is efficiently converted into stars, leaving an extended 
atmosphere of high angular momentum gas. 

There is good agreement between the TREESPH and 
GRAPESPH runs plotted in Fig. ^ with one exception, 
namely the final specific angular momentum of the stellar 
components of the halo 3 runs with Zcooi ~ 1 (see Sec- 
tion 5.2.1, for a detailed discussion). The general conclusion 
from both the TREESPH and GRAPESPH simulations is 
the same however: the stellar systems that form in simu- 
lations with high redshift cooling experience an 'angular- 
momentum catastrophe' and end up with specific angular 
momenta that are between one or two orders of magnitude 
smaller than the specific angular momenta of real disc galax- 
ies. In contrast, if cooling is suppressed until Zcooi ~ 1, stel- 
lar discs can form with angular momenta that are within 
the range observed for real galaxies. For example, compar- 
ing the Zcooi = 1 results plotted in Fig. ^ with the results of 
Figure 1 shows that the stellar systems that form in haloes 
5, 3 (in the TREESPH run) and 2 have similar specific an- 
gular momenta to real disc galaxies. They also have smaller 
specific angular momenta, by a factor of about 4, compared 
to their parent dark matter haloes at z = 0, in agreement 
with the results plotted in Figure 1. These results provide 
powerful evidence that late disc formation is essential if we 
are to explain their angular momenta. 

However, suppressing cooling until ^cooi = 1 does not 
guarantee that disc systems will form by the present day. 
Thus, even with cooling suppressed until such a recent 
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Figure 14. Specific angular momenta at ^ = versus circular speed Vc for stars within 20 kpc (plotted as stars) gas (plotted as triangles) 
and dark matter (plotted as circles) at the virial radius. TREESPH results are shown by the filled symbols and GRAPESPH results by 
the open symbols. For the stellar components, 4— pointed stars show results for runs with Zcool = 4, 5— pointed stars for Zcool = li and 
6— pointed stars for Zcool = 0.6. Crosses are for Halo 2 stars evolved with cooling suppressed by g{z). The abscissa is the circular velocity 
of the halo at the viral radius. The solid line is equation pi 



epoch, the stellar remnants that form in Haloes 1, 4 (and 
Halo 3 in the GRAPESPH run) have much lower specific an- 
gular momenta than real disc galaxies. This is because these 
dark matter haloes evolve considerably between z = 1 and 
the present (nearly doubling their masses within the virial 
radius); star formation in these haloes occurs in sub-units 
which merge between z = 1 and the present day, losing an- 
gular momentum in the process. The formation of discs thus 
depends sensitively on the internal structure and merger his- 
tory of their parent dark matter haloes at recent epochs. 

5.4 Colour Evolution 

The star formation histories for the five models with Zcooi = 
1 are shown in Fig. |l5[ As explained in Section 5.2.2, the 
star formation histories in the TREESPH and GRAPESPH 



runs are slightly different, with more net star formation in 
the TREESPH models. However, in both sets of models, the 
star-formation histories are quite well approximated by the 
form 

M4t') = a(l - exp i-f3t')), (20) 

with values of a and /3 as listed in Fig. ^ In equation |20| 
t' denotes the time measured from z — 1 when cooling is 
switched on. In these models, there is no star formation by 
construction at z > 1, and so the entire present day galaxy 
forms since z = 1. The star formation rates decline once star 
formation begins. In the extreme cases shown in Fig. the 
star formation rates vary between M « YOMq/jt at z = 1 
to M ~ 12M0/yr at 2 = for Halo 1 (which does not form a 
final disc) and M « 28Mq /yr at 2: = 1 and M « 1.5Mq /yr 
at 2; = for Halo 5 (which does form a stellar disc). 
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In the previous subsection, we have argued that high 
angular momenta of disc galaxies requires that they formed 
at late epochs. Thus, the question arises as to whether such 
extreme evolution as implied by Fig. |l^ could be compati- 
ble with observations of disc galaxies at redshifts z ~ 0.5-1 
(Lilly et al. 1997, Vogt et al. 1997a, b). We have therefore 
computed the rest frame absolute magnitudes and colours 
for Haloes 1 and 5 (spanning the range of star formation 
histories plotted in Fig. using the Bruzual and Chariot 
(1993) population synthesis code. The results are plotted 
in Fig. |l6| for a Salpeter (1955) initial stellar mass func- 
tion extending from O.IMq to 125M0. These models show 
that there is weak evolution in the K-band, once star for- 
mation is underway. The rest-frame B-band absolute mag- 
nitudes brighten by more than a magnitude between z = Q 
and z = 0.8, despite the fact that the stellar disc mass de- 
creases by an order of magnitude between these epochs. This 
is not surprising because the luminosity, especially in the 
B-band, is sensitive to the net star-formation rate at early 
times rather than the underlying disc mass. 

An analysis of recent observations of disc evolution 
(Mao, Mo and White 1998), suggests that at fixed circular 
speed, discs show very little evolution in B-band luminos- 
ity. Comparison with Fig. |l^ suggests that the star forma- 
tion rates in our simulations are probably too high at early 
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Figure 16. Bruzual and Chariot models of the rest-frame ab- 
solute magnitude evolution in the B and K bands, and B — I 
colours, for the star formation rates in Haloes 1 and 5 as plotted 
in Figure |l5| The pairs of lines show the evolution determined 
for the star formation rates from the TREESPH and GRAPE- 
SPH simulations. A standard Salpeter (1955) stellar initial mass 
function was assumed. 



times. In the simple cooling scheme adopted in our models, 
a reservoir of high density gas builds up in the central re- 
gions of the haloes that cools within a dynamical time and 
forms stars once cooling is switched on. It is therefore likely 
that the star formation rates will be less strongly peaked to 
high redshifts in more realistic feedback models and hence 
that the luminosity evolution will be even weaker than that 
shown in Fig. ^ 

The very strong evolution of disc masses with redshift 
should manifest itself as a strong evolution of the disc struc- 
tural parameters (scale lengths) with redshift. Discs at high 
redshift should have smaller scale lengths than at the present 
day. However, the gravitational softenings in our simulations 
are comparable to the disc scale lengths and so we cannot 
reliably predict either their size or how they will evolve with 
redshift. (In most of our runs, the stellar disc half-mass radii 
are about equal to the stellar softening and stay roughly con- 
stant as the discs evolve.) Nevertheless, Mao et al. (1998) 
find evidence that at fixed circular speed, disc scale lengths 
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do decrease as (1 + z)~^. This is compatible both with sim- 
ple analytical models of disc evolution {e.g. Mo et al. 1998) 
and with the conclusions in favour of late disc formation 
presented in this paper. 



6 CONCLUSIONS 

The purpose of this study has been to investigate whether 
disc galaxies can form in a CDM universe with angular mo- 
menta and sizes similar to those of real disc systems. We have 
computed the specific angular momenta of CDM haloes from 
an N-body simulation and compared them with estimates 
for a large sample of disc galaxies. This comparison shows 
that disc galaxies have very nearly the same specific angu- 
lar momenta as those within the virialized region of CDM 
haloes at z « 1. We therefore conclude that in a CDM- like 
model, the gas that formed the disc component must have 
collapsed at recent epochs, very nearly conserving its angu- 
lar momentum during collapse. This picture is capable of 
explaining the present day sizes and distribution of surface 
brightnesses of disc galaxies (see also Dalcanton et al. 1997, 
and Mo et al. 1998). 

However, numerical simulations of disc formation have 
shown that the gas experiences an 'angular momentum 
catastrophe' if it is allowed to collapse early {e.g. Navarro 
& Benz, 1991; Navarro, Fronk & White, 1995a; Navarro & 
Steinmetz 1997). In the absence of feedback processes, the 
gas collapses into subunits at high redshift, which subse- 
quently merge losing most of their angular momentum to 
the dark halo material. The specific angular momenta of the 
resulting gaseous discs are typically two orders of magnitude 
smaller than those of real disc galaxies. 

In this paper, we have tested whether the 'angular mo- 
mentum catastrophe' can be avoided if feedback processes 
prevent the gas from collapsing until late epochs when haloes 
are reasonably smooth and slowly evolving. If this can be 
achieved, it may be possible for the gas to conserve most of 
its angular momentum during disc formation, as assumed in 
the simple analytical models of Fall and Efstathiou (1980) 
and others. Wo have not attempted to model feedback in 
any detail in this paper, but have instead suppressed cool- 
ing until a specified redshift. We have selected five haloes 
with present day circular speeds in the range 150 ^ ~ 270 
km/s from a dissipationless N-body simulation of an n = 1 
CDM model. We have generated higher resolution multi- 
mass initial conditions of these haloes which we have evolved 
with two SPH codes, TREESPH and GRAPESPH, that in- 
clude star formation. 

The results from the two codes arc generally in good 
agreement. Most of the differences are caused by the larger 
softening length of the stellar component of the GRAPESPH 
code which affects the internal structure of the collapsed ob- 
jects and their susceptibility to disruption by mergers. In one 
pair of simulations (runs 11 and 12) the final stellar system 
has a different morphological appearance in the two codes, 
appearing disc-like in TREESPH and spheroidal in GRAPE- 
SPH. However, in all other simulations, the morphologies of 
the final stellar systems are qualitatively similar in the two 
codes. 

All of the runs with cooling suppressed until «cooi = 4 
produced spheroidal systems with final specific angular mo- 



menta of between one and two orders of magnitude lower 
than their dark matter haloes. These systems experience 
the 'angular momentum catastrophe' seen in previous simu- 
lations. However, in the simulations with cooling suppressed 
until Zcooi = 1, stellar discs can form with specific angular 
momenta similar to those of real disc galax;ies. This occurs 
in three out of five haloes simulated with TREESPH and 
two out of five simulated with GRAPESPH. Suppressing 
disc formation until late epochs can therefore solve the 'an- 
gular momentum catastrophe' in some cases, but docs not 
guarantee that a disc will survive to the present day. The 
morphology of the final stellar object in the simulations de- 
pends sensitively on the merging history of the parent halo. 
All of the haloes that wc have simulated have been chosen so 
that they do not merge with a comparable mass system at 
z ^ 1. Nevertheless, Haloes 1 and 4, approximately double 
their mass between z = 1 and z = through a succession 
of mergers and do not maintain stellar discs to the present 
day. Haloes 2 and 5 grow least between z = 1 and z = 
(see Table 1) and these produce the most convincing stellar 
discs in our simulations. 

Our simulations have demonstrated that angular mo- 
mentum evolution during galaxy formation is extremely sen- 
sitive to the thermal history of the gas. This suggests that 
feedback processes can solve the 'angular momentum catas- 
trophe' and are a necessary ingredient in disc galaxy forma- 
tion. The structure of dark matter haloes is also important 
in determining the angular momentum evolution. In fact, a 
number of authors (Katz & Gunn, 1991, Vedel et al. , 1994, 
Steinmetz & MuUer, 1995, Contardo, Steinmetz and Fritze- 
von Alvensleben 1998) have shown that gas approximately 
conserves its angular momentum if it collapses within nearly 
uniform dark matter haloes that grow by monolithic collapse 
rather than hierarchical merging. This is not a realistic solu- 
tion, however, within the context of CDM-like theories. The 
detailed evolution of dark matter haloes in CDM models will 
depend on the cosmology and spectrum of fluctuations and 
these may be important in determining the fraction of disc 
galax;ies at the present day. To address this issue in an unbi- 
ased fashion would require a random choice of haloes to be 
resimulated, unlike that adopted here. 

Our simulations suggest that the disc components of 
galaxies formed at recent epochs, perhaps as late as z ^ 1. 
This is in qualitative agreement with recent observations 
showing that discs are evolving significantly over this red- 
shift range (Mao et al. 1998). A more detailed comparison 
with observations will, however, require higher resolution 
simulations and more realistic models of feedback and star 
formation. 



ACKNOWLEDGMENTS 

We thank Yong-Ik Byun and Ken Freeman for allowing us to 
use their disc photometry parameters prior to publication. 
MLW acknowledges funding by a PPARC postdoctoral posi- 
tion at the University of Oxford. VRE acknowledges Douglas 
Heggie for looking after the GRAPE in Edinburgh and the 
support of a PPARC postdoctoral fellowship. GPE thanks 
PPARC for the award of a Senior Fellowship. 



© 0000 RAS, MNRAS 000, 000-000 



Disc galaxy formation 19 



ApJ, 



REFERENCES 

Allen R.J., Shu F.H., 1979, ApJ, 227, 67 
Barnes J., Efstathiou G., 1987, ApJ, 319, 575 
Benz W., Bowers R.L., Cameron A.G.W., Press W.H., 1990, 
348, 647 

Black J.H., 1981, MNRAS, 197, 553 
Bruzual G., Chariot S., 1993, ApJ, 415, 580 
Byun Y.I., Freeman K.C., 1995, ApJ, 448, 563 
Gen R., 1992, ApJS, 78, 341 

Gontardo G., Stein mctz M.. Fritzc-v on Alvensleben U., 1998, 



ApJ, submitted (isto-ph/9801278) 
Dalcanton J.J., Spergel S.N., Summers F.J., 1997, ApJ, 482, 669 
Davis M., Efstathiou G., Frenk G.S., White S.D.M., 1985, ApJ, 

292, 371 

1976, Nature, 263, 573 
1995, MNRAS, 272, L25 

Davis M., Frcnk G.S., White S.D.M., 1985, ApJS, 



Disney M.J 
Efstathiou G. 
Efstathiou G. 

57, 241 
Efstathiou G. 
Efstathiou G 



Ellis R.S., Peterson B.A., 1988, MNRAS, 232, 431 
Barnes J., 1984, in Formation and evolution of 
galaxies and large structures in the universe. D. Rcidel, Dor- 
drecht, p. 361 

Efstathiou G., Bond J.R., White S.D.M., 1992, MNRAS, 258, LI 
Eke V.R., Cole S., Frenk C.S., 1996, MNRAS, 282, 263 
Eke V. R.. Nava rro J.F., Frenk C.S., 1998, ApJ, in press ( astro- 
ph/ |70807o| ) 

Ewell M.W., 1988, Ph.D. thesis, Princeton University 
Fall S.M., 1983, in Internal Kinematics and Dynamics of Galax- 
ies, ed. E. Athanassoula, (Dordrecht: Reidel), p 391. 
Fall S.M., Efstathiou G., 1980, MNRAS, 193, 189 
Frenk C.S., White S.D.M., Davis M., Efstathiou G., 1988, ApJ, 
327, 507 

Gingold R.A., Monaghan J. J., 1977, MNRAS, 181, 375 

Gunn J.E., 1982, in Astrophysical Cosmology, ed. H.A. Bruck, 

G. Coyne and M.S. Longair, (Vatican Pontificia Acadamia 

Scientiarum), 191. 
Hernquist L., Katz N., 1989, ApJS, 70, 419 



Jenkins A.R. et al. 1998, ApJ, in press (a,stro-ph/9709010) 

Katz N., 1992, ApJ, 391, 502 

Katz N., Gunn J.E., 1991, ApJ, 377, 365 

Katz N., White S.D.M., 1993, ApJ, 412, 455 

Lemson G.. Ka uffmann. G., 1997, MNRAS submitted, ( astro- 

ph/ j710125| ) 
Lilly S.J., et al. , 1997, ( |astro-ph/971206l[ ) 
Lucy L., 1977, AJ, 82, 1013 

Mao S.. Mo H.J.. White S.D.M., 1998, submitted to MNRAS 



(astro-ph/9712167) 



Mathewson D.S., Ford V.L., Buchhorn M., 1992, ApJS, 81, 413 
Mo H.J.. Mao S.. White S.D.M., 1998, submitted to MNRAS, 



(astro-ph/9707344) 



J., Lattanzio J.C., 1985, A&A, 149, 135 
, Benz W., 1991, ApJ, 380, 320 
, White S.D.M., 1993, MNRAS, 265, 271 (NW) 
Frenk C.S., White S.D.M., 1995a, MNRAS, 275, 56 
, Frenk C.S., White S.D.M., 1995b, MNRAS, 275, 



Monaghan J 
Navarro J.F. 
Navarro J.F. 
Navarro J.F. 
Navarro J.F 
720 

Navarro J.F., Frenk C.S., White S.D.M., 1997, ApJ, 490, 493 

Navarro J.F., Steinmetz M., 1997, ApJ, 478, 13 

Okumura S.K., Makino J., Ebisuzaki T., Ito T., Fukushige T 

Sugimoto D., 1993, PAS J, 45, 329 
Peebles P.J.E. 1969, ApJ, 155, 393. 
Phillipps S., Disney M., 1986, MNRAS, 221, 1039 
Quinn T., Katz N., Efstathiou G., 1996, MNRAS, 278, L49 
Salpeter E.E., 1955, ApJ, 121, 161 
Steinmetz M., 1996, MNRAS, 278, 1005 
Steinmetz M., MuUer E., 1995, MNRAS, 276, 549 
Steinmetz M., White S.D.M., 1997, MNRAS, 288, 545 



Sugimoto D., Chikada Y., Makino J., Ito T., Ebisuzaki T., 

Umemura M., 1990, Nat, 345, 33 
TuUy R.B., Fisher J.R., 1977, A&A, 54, 661 

Vedel H., Hellsten U., Sommer-Larsen J., 1994, MNRAS, 271, 743 
Vogt N.P., Forbes D.A., PhiUips A.C., Gronwall C, Faber S.M., 

lUingworth G.D., Koo D., 1997a, ApJ, 465, L15 
Vogt N.P., PhiUips A.C., Faber S.M., Gallego J., Gronwall C., 

Guzman R., lUingworth G.D., Koo D., Lowenthal J.D., 1997b, 

ApJ, 479, L121 
White S.D.M., 1984, ApJ, 286, 38 
White S.D.M., Rees M.J., 1978, MNRAS, 183 341. 
White S.D.M., Efstathiou G., Frenk C.S., 1993, MNRAS, 

1023 

Yepes G., Kates R., Khokhlov A., Klypin A., 1997, MNRAS, 
235 

Zurek W.H., Quinn P.J., Salmon J.K., 1988, ApJ, 330, 519. 



262, 



284, 



© 0000 RAS, MNRAS 000, 000-000 



2 Mpc 




TREESPH Halo 2 GRAPESPH 







. . ^ 


A" 


1 

4'' : * 
^ •*. — 





TREESPH Halo 3 GRAPESPH 




TREESPH Halo 3 



(b) 




.f 








J.-. 





GRAPESPH Halo 3 




TREESPH Halo 3 



GRAPESPH Halo 3 



